rm(list=ls())

# Figure A.3

load("survey_data.rda")
stan.mod.6 <- readRDS("claassen_mod6_mean_estimates.rds")

item.names = levels(as.factor(exp_claassen_input[["data"]][["item"]]))
item.pars = data.frame(
  item_int=round(summary(stan.mod.6, pars="lambda")$summary[,"mean"], 3), 
  item_slope=round(summary(stan.mod.6, pars="gamm")$summary[,"mean"], 3))
rownames(item.pars) = item.names
item.pars[order(item.pars[,2]),]

# adjust intercepts and slopes to account for standardization of theta 
theta.out = rstan::extract(stan.mod.6, pars = c("theta"))[[1]]
theta.raw.mean = mean(as.vector(theta.out))
theta.raw.sd = sd(as.vector(theta.out))
item.pars.std = item.pars
# intercept (a): a^std = a + b*mu
# slope (b): b^std = b*sigma
item.pars.std[,1] = item.pars[,1] + (item.pars[,2] * theta.raw.mean)
item.pars.std[,2] = item.pars[,2] * theta.raw.sd
item.pars.std
arm::invlogit(item.pars.std)

# set up plot
plot.pal = c("#000000", "#000000", "#000000", "#838383", "#838383", "#838383", "#838383", "#838383", "#838383" )
plot.lty = c(1, 5, 4, 1, 5, 4, 4, 4, 4)
names.col = rep(9, length(item.names))
names.col[grep("army", item.names)] = 1
names.col[grep("church", item.names)] = 2
names.col[grep("evdemoc", item.names)] = 3
names.col[grep("strong", item.names)] = 4
names.col[grep("threestate", item.names)] = 5
leg.names = c("Army rule", "Churchill", "Evaluate democracy", "Strong leader", "Three statements", "Other")
proj.names = c("AfroBarometer", "ArabBarometer", "AsiaBarometer", "AsianBarometer", "CSES", "WVS", "EVS","LatinoBarometer",
               "South Asian Barometer","Pew", "Other")

proj.inds = length(names.col)
proj.inds[grep("afrob", item.names)] = 1
proj.inds[grep("arabb", item.names)] = 2 
proj.inds[grep("asiab", item.names)] = 3
proj.inds[grep("asianb", item.names)] = 4
proj.inds[grep("cses", item.names)] = 5
proj.inds[grep("wvs", item.names)] = 6
proj.inds[grep("evs", item.names)] = 7
proj.inds[grep("lb", item.names)] = 8
proj.inds[grep("sasianb", item.names)] = 9
proj.inds[grep("pew", item.names)] = 10
proj.inds = ifelse(is.na(proj.inds),11,proj.inds)

# plot
pdf("fig_a3.pdf", width=6.5, height=5.5)
layout(matrix(1:12, nrow=3, ncol=4, byrow=TRUE), widths=c(1.2,1,1,1), heights=c(1,1,1.2))

# plot 1
par(mar=c(0.1, 3, 2, 0.5), cex=0.7)
curve(arm::invlogit(item.pars.std[proj.inds==1,1][1] + x*item.pars.std[1,2]), -4, 4, type="n", 
      xlab="", ylab="", ylim=c(0,1), yaxs="i", xaxs="i", xaxt="n", yaxt="n", bty="n")
abline(h=c(.25,.5,.75), lty=3, lwd=0.75, col=grey(0.5))
abline(v=c(-2,0,2), lty=3, lwd=0.75, col=grey(0.5))
for(j in 1:dim(item.pars.std[proj.inds==1,])[1]) {
  curve(arm::invlogit(item.pars.std[proj.inds==1,1][j] + x*item.pars.std[proj.inds==1,2][j]), 
        -4, 4, type="l", add=TRUE, lwd=1.75, col=plot.pal[names.col[proj.inds==1]][j], 
        lty=plot.lty[names.col[proj.inds==1]][j])
}
mtext(proj.names[1], side=3, line=0.4, cex=0.75)
axis(2, at=c(.25, .5, .75), labels=c(".25", ".50", ".75"), mgp=c(1.5, .35, 0), las=1, 
     tcl=-0.2, cex.axis=0.78)
box(lty=1, lwd=0.75)

# plots 2-4  
for(i in 2:4) {
  par(mar=c(0.1, 0.8, 2, 0.5), cex=0.7)
  curve(arm::invlogit(item.pars.std[proj.inds==i,1][1] + x*item.pars.std[1,2]), -4, 4, type="n", 
        xlab="", ylab="", ylim=c(0,1), yaxs="i", xaxs="i", xaxt="n", yaxt="n", bty="n")
  abline(h=c(.25,.5,.75), lty=3, lwd=0.75, col=grey(0.5))
  abline(v=c(-2,0,2), lty=3, lwd=0.75, col=grey(0.5))
  for(j in 1:dim(item.pars.std[proj.inds==i,])[1]) {
    curve(arm::invlogit(item.pars.std[proj.inds==i,1][j] + x*item.pars.std[proj.inds==i,2][j]), 
          -4, 4, type="l", add=TRUE, lwd=1.75, col=plot.pal[names.col[proj.inds==i]][j], 
          lty=plot.lty[names.col[proj.inds==i]][j])
  }
  mtext(proj.names[i], side=3, line=0.3, cex=0.75)
  box(lty=1, lwd=0.75)
}

# plot 5
par(mar=c(0.1, 3, 2, 0.5), cex=0.7)
curve(arm::invlogit(item.pars.std[proj.inds==5,1][1] + x*item.pars.std[1,2]), -4, 4, type="n", 
      xlab="", ylab="", ylim=c(0,1), yaxs="i", xaxs="i", xaxt="n", yaxt="n", bty="n")
abline(h=c(.25,.5,.75), lty=3, lwd=0.75, col=grey(0.5))
abline(v=c(-2,0,2), lty=3, lwd=0.75, col=grey(0.5))
for(j in 1:dim(item.pars.std[proj.inds==5,])[1]) {
  curve(arm::invlogit(item.pars.std[proj.inds==5,1][j] + x*item.pars.std[proj.inds==5,2][j]), 
        -4, 4, type="l", add=TRUE, lwd=1.75, col=plot.pal[names.col[proj.inds==5]][j], 
        lty=plot.lty[names.col[proj.inds==5]][j])
}
mtext(proj.names[5], side=3, line=0.3, cex=0.75)
axis(2, at=c(.25, .5, .75), labels=c(".25", ".50", ".75"), mgp=c(1.5, .35, 0), las=1, 
     tcl=-0.2, cex.axis=0.78)
mtext("Proportion supporting democracy", side=2, line=1.5, las=0, cex=0.75)
box(lty=1, lwd=0.75)

# plots 6-8  
for(i in 6:8) {
  par(mar=c(0.1, 0.8, 2, 0.5), cex=0.7)
  curve(arm::invlogit(item.pars.std[proj.inds==i,1][1] + x*item.pars.std[1,2]), -4, 4, type="n", 
        xlab="", ylab="", ylim=c(0,1), yaxs="i", xaxs="i", xaxt="n", yaxt="n", bty="n")
  abline(h=c(.25,.5,.75), lty=3, lwd=0.75, col=grey(0.5))
  abline(v=c(-2,0,2), lty=3, lwd=0.75, col=grey(0.5))
  for(j in 1:dim(item.pars.std[proj.inds==i,])[1]) {
    curve(arm::invlogit(item.pars.std[proj.inds==i,1][j] + x*item.pars.std[proj.inds==i,2][j]), 
          -4, 4, type="l", add=TRUE, lwd=1.75, col=plot.pal[names.col[proj.inds==i]][j], 
          lty=plot.lty[names.col[proj.inds==i]][j])
  }
  mtext(proj.names[i], side=3, line=0.3, cex=0.75)
  box(lty=1, lwd=0.75)
}

# plot 9
par(mar=c(2.5, 3, 2, 0.5), cex=0.7)
curve(arm::invlogit(item.pars.std[proj.inds==9,1][1] + x*item.pars.std[1,2]), -4, 4, type="n", 
      xlab="", ylab="", ylim=c(0,1), yaxs="i", xaxs="i", xaxt="n", yaxt="n", bty="n")
abline(h=c(.25,.5,.75), lty=3, lwd=0.75, col=grey(0.5))
abline(v=c(-2,0,2), lty=3, lwd=0.75, col=grey(0.5))
for(j in 1:dim(item.pars.std[proj.inds==9,])[1]) {
  curve(arm::invlogit(item.pars.std[proj.inds==9,1][j] + x*item.pars.std[proj.inds==9,2][j]), 
        -4, 4, type="l", add=TRUE, lwd=1.75, col=plot.pal[names.col[proj.inds==9]][j], 
        lty=plot.lty[names.col[proj.inds==9]][j])
}
axis(1, at=c(-2, 0, 2),  mgp=c(1, .05, 0), las=1, tcl=-0.2, cex.axis=0.78)
mtext(expression(theta), side=1, line=1.2, cex=0.75)
mtext(proj.names[9], side=3, line=0.3, cex=0.75)
axis(2, at=c(.25, .5, .75), labels=c(".25", ".50", ".75"), mgp=c(1.5, .35, 0), las=1, 
     tcl=-0.2, cex.axis=0.78)
box(lty=1, lwd=0.75)

# plots 10-11  
for(i in 10:11) {
  par(mar=c(2.5, 0.8, 2, 0.5), cex=0.7)
  curve(arm::invlogit(item.pars.std[proj.inds==i,1][1] + x*item.pars.std[1,2]), -4, 4, type="n", 
        xlab="", ylab="", ylim=c(0,1), yaxs="i", xaxs="i", xaxt="n", yaxt="n", bty="n")
  abline(h=c(.25,.5,.75), lty=3, lwd=0.75, col=grey(0.5))
  abline(v=c(-2,0,2), lty=3, lwd=0.75, col=grey(0.5))
  for(j in 1:dim(item.pars.std[proj.inds==i,])[1]) {
    curve(arm::invlogit(item.pars.std[proj.inds==i,1][j] + x*item.pars.std[proj.inds==i,2][j]), 
          -4, 4, type="l", add=TRUE, lwd=1.75, col=plot.pal[names.col[proj.inds==i]][j], 
          lty=plot.lty[names.col[proj.inds==i]][j])
  }
  axis(1, at=c(-2, 0, 2),  mgp=c(1, .05, 0), las=1, tcl=-0.2, cex.axis=0.78)
  mtext(expression(theta), side=1, line=1.2, cex=0.75)
  mtext(proj.names[i], side=3, line=0.3, cex=0.75)
  box(lty=1, lwd=0.75)
}

# plot 12
par(mar=c(2.5, 0.8, 2, 0.5), cex=0.7)
curve(arm::invlogit(item.pars.std[proj.inds==9,1][1] + x*item.pars.std[1,2]), -4, 4, type="n", 
      xlab="", ylab="", ylim=c(0,1), yaxs="i", xaxs="i", xaxt="n", yaxt="n", bty="n")
legend(x=-3.6, y=.8, legend=leg.names, lwd=1.75, col=plot.pal[c(1:5,9)], lty=c(1,5,4,1,5,4), bty="n", 
       cex=0.8, seg.len=3)
box(lty=1, lwd=0.75)

dev.off()
